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ABSTRACT 



The satellite kinematics model in the Naval Wargaming System 
(NWGS) is described and critiqued. Suggestions are made for 
simplifying the existing computer code and for increasing the 
accuracy of the satellite simulation. 



A. Introduction 



The Ml 0 satellite control model (SCM) in the Naval Warfare 
Gaming System (NWGS) serves two primary purposes. First, the model 
performs orbital kinematic computations to position the 
satellite (s) to the current game time. Second, the model selects 
either the radar satellite detection model or the elint detection 
model depending upon the type of satellite being processed. The 
purpose of this report is to critique the orbital kinematics 
section of the SCM and to make recommendations for improvement. 

B. The SCM Kinematics 

Satellite parameter inputs to the SCM, pertinent to the 
kinematics, consist of the satellite orbital period, the time of 
last position update and the latitude, longitude and course at the 
time of last position update. The time interval since the last 
position update, t_int, is the current game time minus the time of 
last update. The great circle distance traveled (measured in arc 
length) since the time of last update is t_int/per iod . Using 
spherical trigonometry, the position and course are updated by 
moving the satellite along the great circle determined by its last 
position and course, and the arc distance traveled in t_int. The 
ground position is then corrected by subtracting the amount of the 
Earth's rotation in t_int from the updated longitude. 

In essence, the SCM assumes all satellites to be in 
unperturbed circular Keplerian orbits about the Earth. 
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C, Critique 



The major limitation of the SCM is that all orbits are assumed 
circular. Even when this assumption is valid, the major flaw in 
the existing SCM is the neglect of the precession of the orbit 
caused by the Earth's oblateness. This precession depends upon the 
height and inclination of the orbit and can be as much as 10 
degrees per day which, in turn, can cause a ground position error 
up to 600 n. mi. per day. Neither orbital height nor inclination 
are input parameters to the SCM. 

An additional error, caused by the neglect of the Earth's 
oblateness, is in the computed latitude of the satellite's ground 
track. At current game time the computed (geocentric) latitude can 
be in error by as much as 11.5 n. mi. from the true (geographic) 
latitude at latitudes near 45 degrees. 

D J _aununa£y_.Q£_Ba£aiPmfindatiQas 

The sophistication of the SCM model can be increased in 
several steps. An overview is given below. Details are given in 
Section E. 

1. The programming of the existing SCM can be simplified by 
using alternate formulas from spherical trigonometry. 

2. By including the satellite semimajor axis (which can be 

derived from its orbital period) and the orbital 

inclination, the existing SCM can be easily modified to 

include the orbital precession caused by the Earth's 
oblateness. Without this inclusion, a ground position error 
of 600 n.mi. per day can occur. 
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3. To include noncircular orbits it will be necessary to 
introduce all of the satellite orbital parameters. The 
computations involve the solution of Kepler's equation, the 
transformation of the in-plane coordinates to an inertial 
system, and the reduction of the coordinates to geocentric 
latitude and longitude. 

4. The satellite's geocentric latitude can be converted to 
geographic latitude to account for the Earth's oblateness. 
This conversion can eliminate an 11.5 n.mi. ground position 
error in the N-S direction. A by-product is the 
determination of the satellite's height. 

The computation of the satellite position using Recommenda- 
tions 3 and 4 are time consuming and could degrade real time 
operation of the existing NWGS. It would, however, be possible to 
precompute and store the satellite positions in a table prior to 
running the NWGS. 

S. Detailed Recommendations 

1. Add a two argument quadrant arctangent function, called 

here 'qatanbf', to the atrigf include file. Details of the 
qatan function are in the appendix. Then, the Ml0_satellite_ 
control_ee module can be modified as follows: 
a. Replace lines 201-202 by 

d_long = qatanbf (sinb„dist * sinb__course , cosb_dist * 
cosb_lat - sinb_dist * cosb_course * sinb__lat) . 
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b. Replace lines 206-220 by 



n_course = qatanbf ( sinb_course * cosb_lat , 

cosb_course * cosb_dist - sinb_dist * sinb_lat) . 
c. With these two modifications, lines 224-231 are no longer 
required. 

These equations take advantage of the general spherical 
triangle [Ref. 1]. Both d_long and n_course, as computed above, 
lie in the range of -1 to +1 bams (-pi to +pi) . The 4-part 
formula [Ref. 2] is used to compute n_course directly. The 
possibility of polar crossing is automatically included. 

2. The formula for the precession of the ascending node of the 

orbit, owing to the Earth's oblateness, is 



N = N 0 - 



(3 J 2 cos i)n(T - T 0 ) 
2a 2 



where N is the longitude of the ascending node at time T, N 0 is 
the longitude of the ascending node at time T 0/ J 2 = 

0.001082616 is the Earth's second harmonic shape parameter, i 
is the orbital inclination, n is the satellite motion and a 
is the semimajor axis. The motion n = 2 /P, where P is the 
orbital period. Further, a and P are related by 



P 2 = 47r 2 a 3 /k 2 , 



where k = 0 . 074366 9161 (e. r .) 3//2 /min. is the Earth's 
gravitational constant. Using these relations, the formula for 
the node may be written in one of several ways: 
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N = Ng - 9.9639(a e /a) 7/2 (T - T 0 ) cos i, 

where a g is the Earth's equatorial radius, N is degrees and T 
is measured in days, or 

N = Ng - 216.74 p" 7/3 (T - T 0 ) cos i, 

where N is in degrees, and P and T are measured in minutes, or, 

N = Ng - 282.33 P~ 7//3 t_int cos i, 

where N is in bams, and P and t_int are in seconds as is the 
custom in the NWGS. 

Using the last equation, line 235 of Ml0_satellite_ 
control_ee can be written as 

sat_long = sat_long + d_long - (bams_sec + prec) * t_int, 

- 7/3 

where prec = 282.83 * P / * cos i. Note that prec is a 

constant for each satellite. 

3. Since noncircular orbits do not give rise to uniform motion 

of the ground track over a non-rotating spherical Earth, the 
method used in M10_satellite_control_ee cannot be used. In 
fact, eccentric orbits can give rise to ground tracks 
containing cusps. 

Whatever procedure is used for the orbital computation, it 
will be necessary to obtain all of the satellite parameters. 
These are: the semimajor axis a, the eccentricity e, the 
longitude of the ascending node Ng, the argument of perigee 



- 5 - 



w 0 , the angle of inclination i, and the mean anomaly M 0 for 
some epoch T 0 , or their equivalent. Because of the complexity 
of the computations, it might be advisable to obtain the 
orbital elements distributed by NORAD and to use simplest of 
their existing programs, the SGP model [Ref. 3]. The NORAD 
elements differ slightly from those given above but they also 
compensate for orbital decay due to atmospheric drag. Another 
alternative would be to use equations such as those given in 
Escobal [Ref. 4]. An outline of the necessary computations is 
given below. 

Set up the global constants: 

k = 0.0743669161 = gravitational constant, 

J 2 = 0.001082616 = second harmonic of Earth, 
a g = 1 e.r. = Earth's equatorial radius, 

Tg = 0.0043752695 rad/min = Earth's rotation rate. 

Input a, e, i, w 0 , N 0 , M 0 , and T 0 . 

Compute the satellite constants: 

p = a(l - e 2 ) , 

n = ka 2/ ^ 2 [l + 1.5 J 2 (l " e ^) (1 ” 1.5 sin 2 i)/p 2 ] , 

w' = 1.5 n(2 - 2.5 sin 2 i)/p 2 , 

N' = - (1.5 J 2 n cos i)/p 2 . 

The time interval, measured in minutes from T 0 , is T - T 0 , 
where T is the current time. Compute: 



- 6 - 



w = w 0 + w* (T - T 0 ) , 

N = N + N' (T - T 0 ) , 

M = mod [n (T - T 0 ) , 2 tt] . 

Solve Kepler's equation for E: 

E - e sin E = M. 

Compute the in-plane position components, X w and Y , and 
the radius vector R: 

X = a (cos E - e) , 
w 

Y w = a (1 - e 2 ) 1/2 sin E, 

R = a ( 1 - e cos E) . 

Reduce the in-plane coordinates to equatorial coordinates: 



P x 
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COS w 


cos N - 


sin 


w 


sin 
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cos 


if 
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cos w 
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sin 
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Compute the right ascension, RA, and declination, DEC, of 



the satellite: 

RA = qatan (Y, X) , 

DEC = arccos(Z/R). 

The geocentric latitude of the satellite = DEC. 



Compute 


the unnormalized longitude: 


K 


= RA - Tg _ T’t, 



where t is the total number of minutes elapsed since the 



year, month 


and day of T„ , and T is the Greenwich sidereal 
1 0 g 


time at the 


initial time Tg . T , measured in degrees, is given 


by 




T 

g 


= 99.6909833 + (36000.7689 + 0.00038708T )T , 

u' u 


where is 


given by 


T 

u 


= (JD - 2415020.0)/36525, 



and JD is the Julian Day Number, 

JD = 367Y - int { 7 [ Y+int ( (M + 9)/12)]/4} 

+ int ( 27 5M/ 9) + D + 1721013.5 + UT/24 , 

where Y is the year (1901 1 Y 1 2099), M is the month (1 1 M l 
12), D is the day (1 £ D £ 31), and UT is the universal time in 



hours. Then 


compute the east longitude L, where 


L 

n 


= mod (I/, 2 tt ) . 
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4. 



A conversion of geocentric latitude, L^, to geographic 
latitude for an oblate Earth can be found in Ref. 4. A slightly- 
improved version is given below. Let 

a g = 6378.135 n. mi. = Earth's equatorial radius, 

f = 1/298.26 = Earth's flattening factor, 

r = a R. 
e 

Let L' = DEL. Then compute 

r c = a e {[l - ( 2 f - f 2 )]/ll - ( 2 f - f 2 ) cos 2 L ' ] } 1/2 ' 

L fc = arctan[tan L'/(l - f ) 2 ] , 

H/r - [1 - U C A) 2 sin 2 (L fc - L'] 1/2 
- ( r c /r) cos (L fc - L ' ) , 

DL' = arcsin[(H/r) sin(L fc - L ' ) ] . 

Next, let L' = DEL - DL ' and recompute from the equation for r Q 
until L' no longer varies (convergence is rapid). L fc is the 
geographic longitude of the satellite and H = r(H/r) is the 
height of the satellite above the oblate Earth. 
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Ap pe ndix 



The quadrant arctangent function, designated qatan, is a 
two argument function of rectangular coordinates x and y 
which returns an arctangent value in the range of -it to +rr 
depending upon the quadrant of x and y. In FORTRAN this 
function is called AT AN 2 . It can be computed as follows: 

^ c 

qatan (x,y) = atan{y/[x + e*t ( x— 0 ) ] } 

+ iT*t(x < 0)*[sign(y) + t(y £ 0)], 

where 

e = smallest number representable on the computer, 



t(z) = 1 


when z 


is 


true, 


t(z) = 0 


when z 


is 


false , 






+1 


if 


y > 3, 


sign (y) 


= 


0 


if 


y = 0 , 






-1 


if 


y < 0 . 
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